A density functional theory study of the confined soft ellipsoid fluid 
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A system of soft ellipsoid molecules confined between two planar walls is studied using classical 
Density Functional Theory (DFT) . Both the isotropic and nematic phases are considered. The excess 
free energy is evaluated using two different Ansatze and the intermolecular interaction is incorporated 
using two different direct correlation functions (DCF). The first is a numerical DCF obtained from 
simulations of bulk soft ellipsoid fluids and the second is taken from Parsons-Lee theory. In both 
the isotropic and nematic phases the numerical DCF gives density and order parameter profiles in 
reasonable agreement with simulation. The Parsons-Lee DCF also gives reasonable agreement in 
the isotropic phase but poor agreement in the nematic phase. 



I. INTRODUCTION 

The behaviour of fluids near surfaces and interfaces 
has attracted great interest in recent years. The pres- 
ence of a surface breaks the translational symmetry of 
the fluid, leading to behaviour radically different to that 
in the bulk. This is particularly true for liquid crystalsi 
Many applications for liquid crystalline materials rely on 
the ability to manipulate the preferred alignment through 
the action of an external field. This behaviour is strongly 
influenced by coupling between the LC and bounding sur- 
faces. 

As may be expected from above, there have been 
many studies of liquid crystal systems near surfaces. 
Experimental studies have been performed using meth- 
ods such as scanning tunnelling microscopy 2 , atom force 
microscopj^ and NMR4 among others. There have also 
been simulation studies, using lattice models^, hard&i or 
soft single site models^, and a few studies using atomistic 
models^. It has also been studied by the all commonly 
used theoretical methods in liquid crystal science. Both 
elastic theory and Landau-de Gennes^S theory have been 
applied to this. However, the former makes some as- 
sumptions (e.g. slow director variation) that are invalid 
near a solid substrate, while the latter introduces phe- 
nomenological parameters not that are not easily related 
to microscopic properties. Lastly there have been several 
studies using density functional theory (DFT)&Liiii2ii^ 
or integral equation methodsi^i^. 

In this paper we study a liquid crystal fluid near a solid 
substrate using density functional theory (DFT)i£»i£. 
DFT in principle provides a route to the thermodynamic 
properties of a fluid from knowledge of its microscopic 
(molecular) properties. The intermolecular potential is 
incorporated through the excess free energy functional. 
While the exact form of this is unknown a number of ap- 
proximations to this are commonly used, such as Onsager 
theory^ and Parsons- Lee theorjii2i22i2i. 

Previous DFT calculations for liquid crystalline sys- 
tems have been stymied as the exact direct correlation 
function (DCF) is unknown. Commonly this is approxi- 
mated by the Mayer /-function 

/(ri2,ui,u 2 ) = exp{-/3F(ri2,Ui,u 2 )} - 1 (1) 

where V(ri2, Ui, u 2 ) is the intermolecular potential and 

a — 1 ia the 



tial, for which DFT calculations have mostly been limited 
to, the integral of the Mayer function gives the excluded 
volume between two particles. 

However, recently the DCF has been calculated for 
the soft ellipsoid model for a number of state points in 
both the nematic and isotropic phase22*2iS4. It has been 
shown to reproduce the bulk properties of the nematic 
fluid such as the elastic constants 2 ^,. Thus it is hoped 
that this may provide a good description of the structure 
of the inhomogeneous fluid. 

In this paper we consider one of the simplest situa- 
tions: a fluid of ellipsoidal molecules interacting with a 
structureless wall. Both isotropic and nematic phases 
are considered. The paper is arranged as follows: in the 
next two sections the theory (Sec. 2) and computational 
method are outlined (Sec. 3). The results are then given 
in Sec. 4 and a short summary follows in Sec. 5. 



II. THEORY 
A. Density Functional Theory 

For a system of uniaxial molecules the Grand Potential 
can be written as 

pn [p(r, u)] = (iF ld [p(r, u)] + (3F ex [p(r, u)] (2) 
+ (3 J drdu V ext (r, u)p(r, u) - [3p J drdu p(r, u) 

where p(r, u) is the orientationally dependent single par- 
ticle density distribution, V ext {r,u) is the external po- 
tential, and /i is the chemical potential. Fid [p( r , u )] and 
F ex [p(r, u)] are the ideal and excess free energies, respec- 
tively. The ideal free energy is given by 

(3F ld [p(r, u)] = J drdu p(r, u) [log (A 3 p(r, u)) - l] , 

(3) 

where A is the thermal de Broglie wavelength. The ex- 
cess free energy is in general unknown. Here we use two 
Ansatze for F ex [p(r, u)] . In the first F ex [p(r, u)] is in the 
spirit of Onsager theorjiiLi^ 

PF ex [p(r,u)}= (4) 
-i [ dr 1 du 1 dr 2 du 2 c(r i2 , ui, u 2 Wri, ui)p(r 2 , u 2 ), 
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where c(ri2, 111,112) is the direct correlation function. 

In the second F ex [p(r, u)] is taken to be a density ex- 
pansion around the homogenous fluid with density po( u ) 
truncated at the 2nd order terrn^ 

(3F ex [p(r,u)] = - J driduic (1) (ui) [p(r 1; Ui) - p {u)] 

~~ \ J rfr i rfu i dr 2 du 2 c(ri2,Ui,U 2 ) (5) 
x {p(ri,ui)-/?o(ui)}{/3(r 2) u 2 ) -p (u 2 )}, 

where po( u ) is the density for the bulk, homogenous fluid 
and c^'(u) is the first order direct correlation function. 
This can be identified with the excess chemical potential 
( 0p ex = -cW)S. 

For fixed external and chemical potentials, the equilib- 
rium single particle density is that which minimises the 
Grand Potential. This is a solution to the Euler-Lagrange 
equation 

SQ [p(r,u)] 



£p(r,u) 







(6) 



log {A 3 p(r, u)}-P(i(u)+I3V(t, u)+ 



S{(3F ex [ P (r,u)]| 
5p{r, u) 



= 

The chemical potential is found from the Euler- 
Lagrange equation of the bulk, homogenous fluid, i.e. 
substituting p(r, u) = po(u) into Eq. [7|for both Ansatze 
for the excess free energy. 

The external potential is given by a repulsive Lennard- 
Jones potential 



V ext {z) 



z+0.5 



z+0.5 
, otherwise 



,if z < 2t 



0.5 



(8) 



This potential acts upon the molecular centres, thus giv- 
ing rise to homeotropic (normal) alignment at the sur- 
face. 



B. Direct Correlation Function 

The direct correlation function (DCF), c(ry, Uj, Uj), is 
the central quantity in many theories of liquids and liq- 
uid crystals. For liquid crystalline systems the simplest 
approximation to the DCF is that of Onsager where the 
DCF is taken to be the Mayer /-function. This corre- 
sponds to truncating the virial expansion at the second 
term, and thus is valid only for low densities and large 
elongations. 

In this paper, two different forms of the DCF are con- 
sidered. The first is a numerical DCF calculated from 
simulations of soft ellipsoid molecules^Si^. In principle 
this should be an exact representation of the interactions 
within the fluid. 

The second form is taken from Parsons-Lee 
theorj*i2i2LSi. This corresponds to a resummed 
virial series truncated, as in Onsager theory at the 
second order term. Thus the Parsons-Lee DCF is the 
Mayer /-function multiplied be a density dependent 
prefactor 



where rj = pv mo i is the packing fraction. To calculate the 
Mayer function the intermolecular potential is required. 
For the soft ellipsoids considered in this paper, the po- 
tential is given by 



V(ri2,ui,u 2 ) = 4 



00 



?'i2 - er(fi2,Ui,u 2 ) + <t 





00 



rij - <7(fi2, Ui, U 2 ) + (Jo 



where 

cr(fy , Ui, Uj) = (To 

X f (f y-Ui 



, (10) 



(11) 



1 + XUj.Uj 



1 - XUi-Uj 



and x — (k 2 ~ 1 )/ (^ 2 + 1 ) > where n is molecular length to 
breadth ratio, equal to 3 for the system considered here. 

In both theories the DCF is taken to be that of the 
bulk, homogenous fluid (i.e. the density in the prefactor 
in Eq. 9 is fixed to be the bulk density). Extension 
of Parsons-Lee theory to inhomogenous densities (i.e. a 
prefactor that depends on the spatially varying density) 
is possible although non-triviai2£i2L2£ and application of 
this is beyond the scope of this paper. 



III. METHOD 

Following AllenSiSS, angular dependent functions are 
expended in spherical harmonics. The single particle den- 
sity and its logarithm become 



P(r, u) 



l,m 



(12) 



and 



log{A 3 p(r,u)} =J2plm(r)Y lm (u). (13) 
Similarly the direct correlation function is expanded 

as 

c(ri2,ui,u 2 ) = (14) 
ci 17ni i 2m2 i m (ri2)Yi imi (ui)Y; 2m2 (u 2 )Yi m (f 12 ). 

In the above sums, < I < l m ax and —l<m<l. 
Symmetry dictates that mi + r«2 + m = and that Zi, 
I2, and I are even. 

For a homogenous fluid confined between two struc- 
tureless walls with the wall normal in the z direction the 
density and its logarithm are functions of z only. 

For the first Ansatz to the excess free energy, substitut- 
ing Eq. ^1 Eq. ^1 and Eq. ^]into the grand potential 
Eq. and performing angular and x and y integrations 

onves 
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pn [p(z,u)} 



J d,Z ^2p lm (z) \jpi m {z) - V47T^o| +y^/3 J 

l,m l,m 

~ ^2 9 / dZl I dZ2 C hhm(\zi ~ Z 2 \)pi im {z 1 )pi 2 rn{z2)- 
17 ^ JO J 



dz Vlm(z)plm(z) ~ V47T/3/1 / dz pim(z) 



(15) 



I 

The Vi m (z) are the components of the spherical har- 

monies expansion of the external potential ( V(z,u) = C hhm (z) = W(2/ + 1)tt fir rc hml2mlQ {r)Pi (-) (16) 
Hi,m^ r im{z)Yi m (\i)). For the potential given above, the ^ J z \rJ 

only nonzero term is Vqq(z) = \f^KV ex t(z). 

C hhm (z) is the integral of the DCF over r and <j> for wh ere P t {x) is the Zth Legendre polynomial, 
a fixed z, given by Similarly the Grand Potential for the second Ansatz 

becomes 



/ dz ^p lm lpi m (z)- V4tt5io\ +Y]f3 dz Vi m (z)pi m (z) - V] (3pi m dz p im (z) 
Jo i } Jo , Jo 

u Z,m /,m t,m 

_ /J / / rfz 2 C; l/2m (|zi - Z 2 |)p; im (zi) [pi 2 m(^) - 2^ m ] 

U^.m- 70 • / -°° 

^1 / dz 2 C h l 2m p hm pi 2 rn (17) 

Z ; I ~~ J J-oc 



where pi m = fi8 m o and /; are the spherical harmonics 
coefficients of the logarithm of the bulk ODF. Note that 
here p includes contriubtions from the linear term in the 
excess free energy Eq. [SJand thus depends on orientation. 

The equilibrium particle densities are found by tab- 
ulating the grand potential on a regular grid and then 
numerically minimising Eq. 1151 or Eq. 1171 with respect 
to the pi m (z). This numerical minimisation is performed 
using the conjugate gradients method^. The calcula- 
tions were taken to be converged when the fractional 
change in energy between iterations was less than 0.001 
(i.e. E new — E old < 0.001). Spatial integrals where evalu- 
ated using the trapezium rule with a grid spacing of 0.04, 
while angular integrations were performed using Gauss- 
Legendre quadrature with 800 to 6400 points. 

Once the equilibrium orientationally dependent den- 
sity p{z, u) is known, the bulk density can be calculated 
by integration over the angular variables, i.e. 



J du p(z,u) 
Attp oq (z). 



Po(z) 

(18) 

Similarly the 2nd rank order parameter profile, P2(z), can 
be calculated fromi^ 

/ dup(z, u)P 2 (u) 



Pi{z) 

P2(z) 



J dup(z,u) 

1 P2o{z) 



(19) 



IV. RESULTS 
A. Density and Order Parameter Profiles 

1. Isotropic Phase 

First we present results for calculations performed in 
the isotropic phase, with p* = 0.24 and T* = 0.5. Shown 
in Fig. ^ are the density and order parameter profiles 
obtained from DFT calculations. Shown for comparison 
are profiles found from molecular dynamics simulations 
of a confined system of soft ellipsoids. 

As can be seen from Fig. QJi the density profiles found 
from DFT calculations using Eq. 0] and Eq. [5] are virtu- 
ally indistinguishable. The DFT profiles share the same 
gross features as those obtained from simulation. There 
is a strong maximum at about z = 0.68, followed by a 
minimum at about z = 1.44. There is also a weak max- 
imum at about z — 3.4 in the DFT profiles and z = 3 
in the simulation data. Both the DFT and simulation 
profiles rapidly tend towards the bulk density. 

The most obvious difference between the DFT and sim- 
ulation results is the height of the first peak. In the pro- 
file obtained from the DFT calculations the height of the 
first peak is about 2.02, while in the simulation data the 
it is approximately 1.61. 

Shown in Fig. are the order parameter profiles ob- 
tained from DFT calculations and simulation. As can 
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FIG. 1: (a) Density profiles for isotropic fluid with p* — 
0.24 and T* — 0.5. Solid line is the simulation data, dashed 
line shows the data obtained from DFT calculations using 
Eq. |l|and the numerical DCF, and the dotted line shows the 
data obtained from DFT calculations using Eq. and the 
numerical DCF. 

(b) Order parameter profiles for isotropic fluid with p* — 0.24 
and T* = 0.5. Same symbols as in (a). 



be seen the agreement between the two profiles here is 
not as good as for the density. Both the DFT and sim- 
ulation data have peaks at about z — 0.64. This peak 
arises due to the effect of the surface on the ordering of 
the molecules. The heights of these peaks are both about 
P2 = 0.64. This initial peak is then followed by a decrease 
in the order parameter. 

In the DFT calculations the order parameter then be- 
comes negative with a minimum at about z — 2.4. The 
profile found from simulation has a minimum at about 
the same distance, although in this case the order pa- 
rameter remains positive. Another maximum is present 
at about z = 3.3 in the simulation data, this is also 
present in the profiles obtained from DFT calculations, 
although of much smaller magnitude. Both profiles then 
rapidly decay to a constant value close to indicating 
isotropic ordering in the bulk of the cell. 

Shown in Fig. [2 are the density and order param- 
eter profiles obtained from DFT calculations using the 
Parson-Lee DCF. As for those obtained using the nu- 
merical DCF, the density profiles (Fig. using Eqs. 
01 and [S] are essentially identical and in good agreement 
with the profile obtained from simulation. Again a major 
discrepancy between the DFT and simulation profiles is 
the height of the first neak. Tn this case however, the 
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FIG. 2: (a) Density profiles for isotropic fluid with p* == 
0.24 and T* = 0.5. Solid line is the simulation data, dashed 
line is data obtained from DFT calculations using Eq. gland 
the Parson-Lee DCF, and the dotted line shows the density 
profile obtained from DFT calculations using Eq. |^]and the 
PL DCF. 

(b) Order parameter profiles for isotropic fluid with p* — 0.24 
and T* = 0.5. same symbols as in (a). 

peak height is underestimated by the DFT calculations, 
1.23 compared to 1.61 for the simulation profile. 

The order parameter profiles are shown in Fig. [21 As 
can be seen the first peak in the order parameter profiles 
obtained from DFT calculations are much smaller than in 
the simulation profile and those found from DFT calcula- 
tions using the numerical DCF. For molecules with length 
to width ration of 3, Parsons-Lee theory predicts an 
isotropic-nematic transition at about p — 0.32 21 , whereas 
for simulation the transition is at p = 0.284. Thus for 
Parsons-Lee theory it is further from the nematic phase, 
and hence at lower order, than in simulation. 



2. Nematic Phase 

Shown in Fig. |3| are the density and order parameter 
profiles for the nematic phase (p* — 0.30 and T* = 0.5) 
calculated using found from DFT calculations and Monte 
Carlo simulation. In the simulations, the cell seize was 
adjusted so that the average density in the cell bulk was 
0.30 to aid comparison between simulation and theory. 
As can be seen the profiles obtained from DFT calcula- 
tions using Eqs. 01 anal are m better agreement with 
each other than with the simulation profiles. 

The densitv nrofiles are shown Fie. DTK. As in the 
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FIG. 3: (a) Density profiles for nematic fluid with p* = 0.30 
and T* = 0.5. Solid line is the simulation data, dashed line 
is the data obtained from DFT calculations using Eq. ^]and 
numerical DCF, and the dotted line shows the data obtained 
from DFT calculations obtained from DFT calculations using 
Eq. |S|and the numerical DCF. 

(b) Order parameter profiles for nematic fluid with p* = 0.30 
and T* = 0.5. Same symbols as in (a). 



FIG. 4: (a) Density profiles for nematic fluid. Solid line is the 
simulation data (p* — 0.30 and T* = 0.5), dashed line is data 
obtained from DFT calculations using Eq. 21 and Parson- Lee 
DCF (p* = 0.35 and T* = 0.5), and the dotted line shows 
the data obtained from DFT calculations using Eq. and 
Parson-Lee DCF (p* = 0.35 and T*0.5). 

(b) Order parameter profiles for nematic fluid. Same symbols 
as in (a). 



isotropic case, there is an initial peak in the density pro- 
file at about z — 0.6. There are then further peaks. In 
the DFT profiles these subsequent peaks are at z = 3.5, 
z = 6.3, and z = 9.4. In the simulation profile, these are 
at z = 3.2, z = 5.6, and z = 7.9. A similar difference has 
been seen in comparison between DFT calculations and 
simulations for simple fluids^i^. In£I this was explained 
due to the imposition of homogeneity in the transverse 
(x and y) directions in the theoretical calculations. 

Aside from the positions of these secondary peaks, the 
other noticeable difference between the DFT and simu- 
lation profiles are in the heights of these peaks. As in 
the isotropic case the height of the initial peak is over- 
estimated compared to simulation. The initial peak in 
the DFT calculations has a height of 3.9, while in the 
simulation profile its height is 2.9. In contrast to the 
isotropic phase, the secondary peaks in the DFT profile 
are more noticeable than those in the simulation profile. 
These secondary peaks are also stronger in the DFT pro- 
file calculated using Eq. 0|than those in the profile calcu- 
lated using Eq. [5] These peaks indicate surface induced 
layering 1 . This has been observed in by X-ray reflectiv- 
ity measurements, both at the nematic-solid interface as 
studied here and at the free nematic interface. 

The second rank order narameter nrofiles are shown 



in Fig. |2t> • Again these show a strong initial peak at 
approximately z — 0.6. This is then follow by a series 
of peaks before decaying to a bulk values as for the den- 
sity profiles. These peaks are at approximately the same 
positions as the peaks in the density profiles. As for 
the density profile the secondary peaks are stronger than 
those in the simulation profile. In the centre of the cell 
the DFT order parameters reach constant values. For 
Eq. 0| the order parameter in the centre of the cell is 
0.65, whereas for Eq. EI (-^) = 0.69. These values are 
close to the bulk of 0.69. The order parameter in the 
centre of the cell from simulation is 0.74, higher than the 
bulk value. 

One major difference between the DFT and simulation 
profiles is the decay of the initial peak. In the DFT profile 
there is a rapid decay from the first peak that is not seen 
in the simulation profile. This can also be seen, albeit 
not as dramatically, in the order parameter profiles for 
the isotropic fluid (Figs. QJandEJ. 

Shown in Fig. 0] are the density and order parameter 
profiles calculated using the Parson-Lee DCF and those 
from simulation. For this model, Parson-Lee theory pre- 
dicts that a system with p* = 0.30 is isotropioSi, these 
calculations were performed with a higher bulk density 

( n* = 3f»1 While this difference makes detailed com- 
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parison impossible, it is small enough for qualitative com- 
parison. 

As can be seen the density and order parameter pro- 
files calculated using Eqs. 0]and[5]are essentially identical 
and are similar to the profiles calculated in the isotropic 
phase using the Parsons-Lee DCF. This can be under- 
stood, as the Parsons-Lee DCF is proportional to the 
Mayer function at all densities, so aside from a constant 
the Parsons-Lee DCF is the same in both the isotropic 
and nematic phases. This leads to a lack of structure in 
the nematic density and order parameter profiles com- 
pared to those from simulation and DCF calculations us- 
ing the numerical DCF. 



B. Variable Cell Width 

In this section, the effect of varying the width of the 
cell is examined. A fluid confined between two planar 
walls is often used in the evaluation of the depletion force 
between two colloidal particles^. Here the surface free 
energy 7, given by the excess (over bulk) grand potential 
per unit areaii 

= fm\p(z,u)]-(m\po(u)] ^ 

is calculated for cell widths from 4 to 40 a. 

The surface free energy as a function of cell width is 
shown in Fig. [5j In all cases the it reaches a constant 
value for large cell widths (which is the surface tension). 
For large cell widths the walls become isolated from each 
other with a layer of bulk fluid intervening between the 
surface regions of fluid. At smaller cell widths there is 
no bulk fluid layer so the surface layers interact directly, 
leading to variation in the surface free energy with cell 
width. This behaviour is generally oscillatory where the 
distance between peaks is of the order of the molecular 
length (2-3 a). 

Shown in Fig. [SJi is the surface free energy as a func- 
tion of cell width calculated using the numerical DCF. 
As can be seen the surface free for the nematic fluid is 
of larger magntude than that of the isotropic fluid and 
shows much stronger variation with cell width. For the 
isotropic fluid it reaches a constant value for cell widths 
greater than about 16 a. In contrast in the nematic fluid 
it only approaches a constant value for cell widths greater 
than about 25 a. This results from the larger range of the 
DCF in the nematic phase compared to in the isotropic 
phase. 

Figure \Sjp shows the surface free energy as a function 
of cell width calculated using the Parsons-Lee DCF. As 
for the numerical DCF the surface free energy in the ne- 
matic phase is of greater magnitude and shows stronger 
variation with cell width than in the isotropic phase. In 
both the nematic and isotropic phases the it is constant 
for L greater than about 18 a. This is in contrast to the 
heha.vinur nf the surface tension for the numerical DCF. 
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FIG. 5: (a) Variation of surface free energy with cell width 
calculated using the numerical DCF. The solid line with cir- 
cles shows the energy of the isotropic fluid calculated using 
Eq. the dashed line with circles shows the energy of the 
isotropic fluid calculated using Eq. |KJ the solid line with 
squares shows the energy of the nematic fluid calculated using 
Eq. [I] and the dashed line with squares shows the energy of 
the nematic fluid calculated using Eq. 03 
(b) Variation of the surface free energy with cell width calcu- 
lated using the Parsons-Lee DCF. The solid line with circles 
shows the energy of the isotropic fluid calculated using Eq. 0] 
the dashed line with circles shows the energy of the isotropic 
fluid calculated using Eq. |S] the solid line with squares shows 
the energy of the nematic fluid calculated using Eq. H and 
the dashed line with squares shows the energy of the nematic 
fluid calculated using Eq. |S] 



V. SUMMARY 

In this paper the structure of a fluid of soft ellipsoids 
near a soft wall is determined using density functional 
theory calculations. The calculated density and order 
parameters have been compared to simulation results for 
the same system. Within the DFT calculations the excess 
free energy was obtained using both a numerical direct 
correlation function and a well known approximation to 
the DCF. Two different Ansdtze were used for the excess 
free energy, giving very similar results to each other. 

There is qualitative agreement between the profiles 
found using the DFT method and those from simulation. 
The agreement is better in the isotropic phase, than for 
the nematic. In the isotropic case, the DFT calculations 
give profiles that appear to be less structured than those 
from simulation. This situation is reversed in the nematic 
nha.se. 
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For the isotropic phase, the Parsons-Lee DCF gave pro- 
files similar to those from the numerical DCF. However, 
in the nematic phase the Parsons-Lee DCF gave profiles 
that were much less structured than those from the nu- 
merical DCF and simulation. 

The variation in the surface free energy as a function of 
cell width has also been examined. For larger cell widths 
it is a constant, whereas for small cell widths it reflects 
strong oscillatory interactions between the two surfaces. 

In summary, the simple density functionals used in 
this paper which are based on a second order density 
expansion of the excess free energy, already give profiles 
which are in qualitative agreement with simulations. We 
still observe quantitative deviations, especially in the ne- 



matic phase. Therefore it will be interesting to compare 
these results with profiles obtained from more sophisti- 
cated density functionals, e. g., the Somoza-Tarazona 
functional^ or the recently proposed modified Roscnfcld 
functional^. 
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